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APPENDIX A 



PINECONE.H 



5 #include <stdio.h> 
/include <stdlib.h> 
/include <math.h> 
/include <gl/gl.h> 
/include <gl/device.h> 

10 



EXPRESS MAIL 364370875 US 



WYCaiki 43Sn.pt, W2S/9S (313) 



#define DC_VALUE_RATIO 0.2 

long number_freqs • 4; 

double freq[10] = { 
1.9, 
1.654, 
1.476, 
1.227 

}; 

double phase [10] = { 
2.111, 
0.0765, 
1.32, 
1.0, 
2.38, 
0.73, 
3.0, 
1.1 

}; 
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STAMPIT.C 

/* this program explores the specifics of stamping an image with 'digital 
reticles' for the purpose of 
5 determining registration between a given image and the digimarc signatures, 
i.e., the registration 

requirements for finding the scale and rotation elements of the ubiquitous 
snowy image patterns */ 

10 



/include "pinecone.h" 
/include "freqs.h" 

15 

/define DISPLAY_IT 1 
/define PI 3.141592653589 
/define SQRT2_2 0.707106781186 
/define INITIAL_SCALE 10.0 

20 

unsigned char *img; 
long xdim,ydim; 
Colorindex *disp; 
Colorindex *scale_disp; 
25 Colorindex *pdisp; 

int helplines ■ 2; 
char *help[] » 
{ 

30 "usage: stamp_it filename xdim ydim channels \n n , 

"\n stamp_it modifies the input image to include digital reticles and 
outputs filename. stamped " 
}; 

35 

int load_scale_disp(void) { 
long i,j; 

pdisp = scale_disp+20; 
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for(j»0; j<15; }++){ 

for(i»l;i<20;i++) { 

*pdisp » (unsigned char)255; 
pdisp +» 20; 

5 } 

pdisp += 20; 

} 

pdisp = scale_disp+20+ ( 35*400) ; 
for(j-0;j<15;j++){ 
10 for(i-l;i<20;i++) { 

♦pdisp = (unsigned char) 255; 

pdisp += 20; 

} 

pdisp 20; 

15 } 
return(O) ; 

} 

int draw_scale( double value) { 
20 long i; 

double dtemp; 

value *=10.0; 

25 memcpy (disp, scale_disp, 40000 ) ; 

dtemp * 20.0 * value; 

if (dtemp > 399.0)dtemp=399.0; 

pdisp = &disp[ 15*400 + (int) dtemp]; 

for ( i*0 ; i<20 ; i++ ) { 
30 *pdisp * (unsigned char) 255; 

pdisp += 400; 

> 

rectwr ite ( 0 , 0 , 399 , 49 , disp ) ; 

35 return (0) ; 
} 
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double calculate_change( double distance , double scale , long which ){ 
long i; 

double value=DC VALUE RATIO * scale; 



5 for ( i«0 ; i<number_f reqs ; i++ ) { 

value += scale * sin { freq[i] * distance + phase(i] + 
(double) which * PI ); 
} 

if(value < 0.0)value * 0.0; 

10 

return (value) ; 
} 



15 

main (argc, argv ) 
int argc ; 

char *argv[] ; 

{ 

20 long i, j ,go-l, button, channels, which, count; 

long temp, increment, last_middle/ 

unsigned char tmp, *pimg, *pimgl, *img_out , *img_ r, *img_b; 
char string [80], outfile(80]; 
FILE *inf; 

25 double change, current_scale * INITIAL_SCALE; 

long gid_img, gid_stamped, gid_scale; 

if(argcl=5) { 

for( j«0; j<helplines; j++)fprintf ( stderr, "%s w , help[j]) ; 
30 exit( 1 ) ; 

} 



xdim = atoi(argv[2] ) ; 
ydim » atoi (argv[ 3 ] ) ; 
35 channels » atoi (argv[ 4] ) ; 

if (channels != 1 && channels 1=3) { 

fprintf ( stderr, "stamp_it : channels must equal 1 for B/W or 
3 for color\n" ) ; 

exit( 1 ) ; 
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} 

if (DISPLAY_IT) disp » calloc(xdim*ydim, sizeof (Colorindex) ) ; 
else disp = 031100(40000, sizeof (Colorindex) ) ; 
5 scale_disp - calloc( 40000, sizeof (Colorindex) ). ; 

img - calloc(xdim*ydim, sizeof (unsigned char) ) ; 
img_out * calloc(xdim*ydim, sizeof (unsigned char) ) ; 
if( Idisp ii lscale_disp ji limg !! limg_out ){ 

fprintf( stderr, "stamp_it : can not allocate space\n w ) ; 
10 exit( 1 ) ; 

} 

if (channels == 3 ) { 

img_r » calloc(xdim*ydim, sizeof (unsigned char) ) ; 
img_b = calloc(xdim*ydim, sizeof (unsigned char) ) ; 
if ( iimg_r | | timg_b ) { 

fprintf( stderr, rt stamp_it : can not allocate space\n" 

exit( 1 ) ; 

} 

} 



15 

) ; 

20 



/* read in binary data into array */ 
25 inf = fopen(argv{ 1 ] , w r M ) ; 

if (i inf) { 

fprintf (stderr, "stamp_it: can't open %s\n" , argv( 1] ) ; 
exit(l); 

} 

30 if (channels « 3){ 

fread(img_r, sizeof (unsigned char) ,xdim*ydim, inf ) ; 
fread(img, sizeof (unsigned char) ,xdim*ydim, inf) ; 
fread(imgjt), sizeof (unsigned char) , xdim*ydim, inf ) ; 

} 

35 else fread( img, sizeof (unsigned char) ,xdim*ydira, inf ) ; 

f close (inf ) ; 



/* flip it */ 
pimg » img; 
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10 



pimgl ■ &img[xdim* (ydim-1) ] ; 

for(i»0;i<(ydim/2) ;i++) { 

for ( j =0 ; j <xdim ; j ++ ) { 
tmp * *pimg/ 
*(pimg++) » *pimgl; 
*(pimgl++) * tmp; 

} 

pimgl -= (2*xdim); 

} 



if (DISPLAY_IT) { 

foreground ( ) ; 
pref size(xdim,ydim) ; 
15 gid_img » winopen( "Original" ) ; 

gf lush( ) ; 

pdisp = disp; 
pimg - img; 

20 for(i=0;i<(ydim*xdim) ;i++) < 

*{pdisp++) = { Color index) * (pimg++) ; 

> 

rectwr ite (0,0, xdim- 1 , ydim- 1 , d isp ) ; 

25 pref size(xdim, ydim) ; 

gid_stamped = winopen( "Stamped Image - . • " ) ; 
rectwrite(0,0,xdim-l,ydim-l,disp) ; 

load_scale_disp( ) ; 
30 pref size (400, 50) ; 



gid_scale » winopen( "Rough scaling..."); 



} 



/* main visual feedback loop */ 
last_middle =0; 
while (go) { 
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/* first diagonal */ 
for(i*0;i<(xdim+ydim-l) ;i++) { 
which * 0; 

/* calculate addition or subtraction to image */ 
5 change » calculate_change( (double) i * SQRT2_2 , 

current^scale, which ) ; 

if{ i < xdim ){ 

pimg = &img[ i*xdim] ; 

pimgl * &img_out[ i*xdim] ; 
10 count = i+1; 

} 

else { 

pimg = &img[xdim*ydim - ydim - xdim + i + 1]; 
pimgl = &img_out[ xdim* ydim - xdim - ydim + i + 

15 1J; 

count = xdim + ydim - i - 1; 

} 

for( j«0; j<count; j++) { 

temp = (long)*pimg + { long) { change+0. 5 ) ; 
20 if (temp > 255) temp ■ 255; 

*pimgl » (unsigned char)temp; 
pimg -= (xdim-1); 
pimgl -« (xdim-1); 

} 

25 } 

/* second diagonal */ 
for ( i=0 ; i< ( xdim+ydim-1 ) ; i++ ) { 
which « 1; 

/* calculate addition or subtraction to image */ 
30 change = calculate_change( (double) (xdim - 1 - i) * 

SQRT2_2 , current_scale, which ); 

if( i < xdim ){ 

pimgl = &img_out[xdim - i - 1]; 
count = i+1; 

35 } 

else { 

pimgl = fiimg_out[ ( i-xdim+1) *xdim] ; 
count * xdim+ydim-i-1; 

> 
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for( j*0; j<count; j++) { 

temp « (long)*pimgl + ( long) (change+O. 5) ; 
if (temp > 255) temp = 255; 
♦pimgl » (unsigned char) temp; 
pimgl += (xdim+1) ; 

} 



10 if (DISPLAY_IT) { 

pdisp = disp; 
pimg = img_ out; 
for(i=0;i<(ydim*xdim) ;i++) { 

*(pdisp++) = (Colorindex)*(pimg++) ; 

15 } 

winset (gid_stamped) ; 

rec twr it e (0,0, xdim- 1 , yd im- 1 , d i sp ) ; 

pdisp = disp; 
20 pimg = img; 

for(i=0;i<(ydim*xdim) ;i++) { 

* (pdisp++) = (Color index) * (pimg++) ; 

> 

winset (gid_img) ; 

25 rec twr i t e (0,0, xd im- 1 , yd im- 1 , d isp ) ; 

winset (gid_scale) ; 
draw_scale(current_scale) ; 

} 

30 

button=0; 
while( ibutton) { 

if( getbutton(LEFTMOUSE) ){ 

current_scale * = 1.15; 
35 button » 1; 

last_middle = 0; 

} 

if( getbutton(RIGHTMOUSE) ){ 

current scale *= 0.85; 



APPENDIX A, Page 97 



WYCzdb 4357t.pi 9/2S/W 013) 



-98 



button * 1; 
lastjniddle * 0; 

} 

if( getbutton(MIDDLEMOUSE) ){ 
5 while( getbutton(MIDDLEMOUSE) ); 

if( lastjniddle ){ 
button 3 1; 
go=0; 

} 

10 lastjniddle * 1; 

} 



15 



/* now re-sign output image because of slight changes to master key *j 



/* flip output image */ 
pimg = img_out; 
20 pimgl = fiimgjaut [xdim* (ydim-1) ] ; 

for(i=0;i<(ydim/2) ;i++) { 

for ( j=0; j<xdim; j++) { 
tmp = *pimg; 
*(pimg++) « *pimgl; 
25 *(pimgl++) = tmp; 

} 

pimgl -= (2*xdim); 

> 

30 /* write out signed image */ 

sprintf (outf ile, -%b. stamped- ,argv( 1] ) ; 
inf = fopen(outf ile, "w") ; 
if(linf) { 

fprintf (stderr, w stamp_it: can't open %s\n" ,outf ile) ; 
35 exit(l); 
} 

if (channels « 3){ 

fwrite( img_r, sizeof (unsigned char) ,xdim*ydim, inf) ; 
fwrite(img_out,sizeof (unsigned char) ,xdim*ydim, inf ) ; 
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fwrite(img_b,eizeof (unsigned char ) , xdim*ydim, inf ) 

} 

else fwrite(img_out,sizeof (unsigned char ) ,xdim*ydim, inf ) ; 
fcloae(inf ) ; 

/* free and clean up */ 
if (DISPLAY_IT) gexit(); 
free(img) ; 
free ( img_out ) ; 
free(disp) ; 
free(scale_disp) ; 
if( channels = 3)< 

f ree( img_r) ; 

f ree(img_b) ; 

} 
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REGISTR.C 

/* this program is a companion to stamp_it, wherein it is given a suspect 
image and it attempts to determine where 

the cross-hatch pattern resides within the suspect image, thereby determining 
the scale, roation and offset 
of the suspect image */ 



10 /include "pinecone.h* 1 
/include "freqs.h" 

/define DISPLAY_IT 1 

/define DISP_POW 0.1 
15 /define MOV_AV 21 /* keep this odd please */ 

/define MAGG_THRE S HOLD 1.7 

/define ANGLE_I NCREMENT (PI/360.0) 

/define START_SCALE_ZONE 50 

/define S CALEBS TART 0.5 
20 /define SCALE_STOP 2.0 

/define SCALE_STEP 0.001 

/define MAX_BLOCK_SIZE 5 

/define PI 3.141592653589 

25 unsigned char *img; 

long xdira,ydim; 

Colorindex *disp; 

Colorindex *pdisp; 

long bits, f ftdim, f ft_size; 
30 float *ar,*ai; 

float wr I 10000] ,wi[ 10000] ; 

int helplines =2; 
char *help[] = 
35 { 

"usage: register filename xdim ydim channels \n", 
"\n register looks at the input image for digital reticles and 
displays registered result " 

}; 
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int shift_array( float *array,int dim) { 
int i,j; 

int dim2 « dim/2; 

int offset = dim2*dim + dim2; 

float *pl, *p2, ftmp; 

for ( i*0 ; i<dim2 ; i++ ) { 

pi = £array[i*dim] ; 

p2 « £array[of f set+i*dim] ; 

for( j=0; j<dim2; j++) { 

ftmp = *pl; 

*pl a *p2; 

*p2 » ftmp; 

pl++;p2++; 

} 

} 

offset » dim2*dim; 
for { i=0 ; i<dim2 ; i++ ) { 

pi = fiarray [dim2+i*dim] ; 
p2 = &array[of fset+i*dim] ; 
for( j=0; j<dim2; j++) { 
ftmp = *pl; 
*pl = *p2; 
*p2 = ftmp; 
pl++;p2++; 

} 

} 

return ( 0 ) ; 

} 

int block_f ilter( float *array,int dim){ 
int i,j,k, 1,12; 
float buffer[2] [4096] ; 

for (i=0;i<dim;i++)buf fer[0] [ i] = array[i]; 
for ( i=l ; i< ( dim-1 ) ; i++ ) { 
i2 = i%2; 

buffer [i2] [0] = array [ i*dim] ; 

buf fer[ i2] [dim-1] = array [ i*dim+dim-l ] 
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for(j-l;j<(dim-l);j++){ 

buffer[i2] [ j ]=0.0; 
for(k»-l;k<2;k++) { 

for(l«-l;K2;l++) { 

bufferf i2] [ array [( i+k) *dim + 



(j+DJ; 



} 

> 

buffer[i2] [ j]/=9.0; 

10 } 

12 = (i-l)%2; 

for< j=0; j<dim; j++)array[ (i-l)*dim + j ] = buf f er [ i2 ] [ j ] ; 

} 

12 = (i-l)%2; 

15 for( j=0; j<dim; j++)array[ (i-l)*dim + j] = buf f er [ i2 ] [ j ] ; 

return ( 0 ) ; 



20 



/* assumes fftdim arrays ar and ai */ 
double get_mag( double x, double y) { 
25 long xoff,yoff; 

float *par,*pai; 

double xdist , ydist , cf , revalue , i_value , value=>0 • 0 ; 

xoff * (long)x; 
30 yoff » (long)y; 

xdist * x - (double) xof f ; 
ydlst « y - ( double ) yof f; 

par ~ &ar{yoff * fftdim + xoff]; 
35 pai = &ai[yoff * fftdim + xoff]; 



cf * (1.0 - xdist) * (1.0-ydist); 
revalue » ( double ) * ( par++ ) ; 
i_value » (double) * (pai++) ; 
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/ 

value » cf * sqrt (r_value*r_value + i_value*i_value) ; 
cf « xdist * (1.0-ydist); 
revalue » (double) *par; 
i — value * (double) *pai; 

value +■ cf * sqrt( revalue* revalue + i_value*i_value) ; 
par +* (fftdim-1); 
pai +- (fftdim-1) ; 

cf « (1.0 - xdist) * ydist; 
revalue * (double) *{par++) ; 
i_value » (double) *(pai++); 

value += cf * sqrt (r_value*r_value + i_value*i_value) ; 
cf = xdist * ydist; 
r_value = (double) *par; 
i_value * (double) *pai; 

value +» cf * sqrt (r_value*r_value + i_value*i_value) ; 
return (value) ; 

} 



main( argc, argv ) 
int argc ; 

char * argv [J ; 
{ 

long i , j , k, go , channels , count , center ; 

long dim, angle_int , total_angles , top_candidate ; 

unsigned char tmp, *pimg, *pimgl, *img_out , *img_r , *img_b; 

char string [80], out file [80]; 

FILE *inf; 

long gid_img; 

double 

x , y , cos8 , s inn , angle , mag, magg [ 5000 ] , magg_ma [ 5000 ] , angle_vector [ 10000 ] ; 

double radius , dtmp, dscale , grey_dif f , highest , f rac , highest_scale, scale 
float *par,*pai; 
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if(argcl-S) { 

for(j«0;j<helplines;j++) fprintf ( stderr, "%s" , help[j]) ; 
exit( 1 ) ; 

> 

xdim » atoi(argv[2 ] ) ; 

ydim ■ atoi(argv(3] ) ; 

channels » atoi(argv[4 ] ) ; 

if (channels 1= 1 && channels l«3){ 

fprintf { stderr, "register : channels must equal 1 for B/W or 
3 for color\n" ) ; 

exit( 1 ) ; 

} 

/* find the next power of two equal to or higher than the highest 
input dimension */ 

if (xdim > ydim) dim = xdim; 
else dim = ydim; 
fftdim « 1; go = 1; bits = 0; 
while ( go ){ 

if{ dim > fftdim ){ 
f ftdim*=2; 
bits++; 

} 

else go = 0; 

> 

if (bits > 12) { 

fprintf ( stderr, "recognize : sorry, this particular program 
only accepts 4K images and less\n" ) ; 
exit( 1 ) ; 

} 

fft_size = fftdim * fftdim; 

disp = calloc(f ft_size, sizeof (Colorindex) ) ; 
img = calloc (xdim*ydim, sizeof (unsigned char) ) ; 
ar = calloc (f ft_size, sizeof (float) ) ; 
ai » calloc(f ft_size, sizeof ( float ) ) ; 
if( tdisp |j limg || lar !! lai ){ 

fprintf ( stderr, "register : can not allocate space\n" ) ; 
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exit( 1 ) ; 

} 

if (channels =*= 3) { 

5 img_r * calloc(xdim*ydim, sizeof (unsigned char) ) ; 

img_b * calloc(xdim*ydim, sizeof (unsigned char) ) ; 
if( limg_r j| iimgjD ){ 

fprintf( stderr, -register : can not allocate space\n" 

) ; 

10 exit( 1 ) ; 

} 

> 

/* read in binary data into array */ 
15 inf = fopen(argv[ 1 ] , "r" ) ; 

if(linf) { 

fprintf (stderr, "register: can't open %s\n" , argv[ 1 ] ) ; 
exit(l); 

> 

20 if (channels »= 3){ 

f read( img_r, sizeof (unsigned char) ,xdim*ydim, inf) ; 
fread( img, sizeof (unsigned char) ,xdim*ydim, inf ) ; 
f read ( img_b, sizeof (unsigned char) ,xdim*ydim, inf) ; 

} 

25 else f read ( img, sizeof (unsigned char) , xdim*ydim, inf ) ; 

f close(inf ) ; 

/* flip it */ 
pimg a img; 
30 pimgl » &img[xdim* (ydim-1 ) ) ; 

for(i=0;i<(ydim/2) ;i++) { 

for( j=0; j<xdim; j++) { 
tmp = *pimg; 
*(pimg++) = *pimgl; 
35 *(pimgl++) » tmp; 

} 

pimgl -» (2*xdim); 

} 
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/* copy image buffer into ar */ 

par * ar; 

pimg * img; 

for ( i»0 ; i<ydim ; i++ ) { 

for{ j=0; j<xdim; j++) { 

*(par++) = (f loat)*(pimg++) ; 

} 

par += (f ftdim-xdim) ; 



/* 2d fft, in place */ 
printf ( w \nforward fft... \n"); 
f ft2d(ar / ai,bits,0,wr,wi) ; 
printf ( "done \n.... "); 
shift_array(ar, f ftdim) ; 
shift_array (ai, f ftdim) ; 
center = fftdim/2; 

block_filter(ar,f ftdim) ; 
block_f ilter(ai, f ftdim) ; 



if (DISPLAY_IT) { 

foreground ( ) ; 
pref size ( f ftdim, f ftdim) ; 
gid_img = winopen( H yahh" ) ; 
gf lush{ ) ; 

dtmp * 

(double) ar [ center* ff tdim+center+1 ] * (double) ar ( center* f f tdira+center+1 ] ; 
dtmp += 

( double ) ai ( center * f f tdim+center+1 ] * ( double ) ai ( center * f f tdim+center+1 ] ; 
dscale = pow(dtmp, DISP_POW); 



pdisp = disp; 
par = ar; 
pai ~ ai; 

for(i=0;i<(f ftdim* f ftdim) ;i++) { 
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dtmp = pow( (*par * *par + *pai * *pai) , 

DISP_POW ) ; 

dtmp *= (255.0 / dscale); 
if (dtmp>255.0)dtmp * 255.0; 
*(pdisp++) =» ( Color index ) { dtmp ); 
par++;pai++; 

} 

rectwrite(0,0, f ftdim-1, f ftdim-l,disp) ; 

> 

/* now search for the gross rotation axes */ 

for(angle = 0. 0, angle_int = 0 ; ang 1 e<P I / 2 ; ang 1 e+ = ANGLE_I NCREMENT / 
angle_int++) { 

coss = cos (angle); 

sinn = sin(angle); 

/* fill radial vector */ 

for (i=0;i<(f ftdim/2) ;i++) { 

radius = (double) i; 

x = (double) center - radius * coss; 

y = (double) center - radius * sinn; 

mag » get_mag(x,y) ; 

x .» (double) center + radius * sinn; 
y = (double) center - radius * coss; 
mag += get_mag(x,y) ; 
magg[i] =* mag; 

} 

/* create moving average */ 
magg_ma[MOV_AV/2] = 0.0; 
for(i=0;i<MOV_AV;i++) { 

magg_ma[MOV_AV/2 J += magg[ij; 

} 

magg_ma[MOV_AV/2] /» ( ( double )MOV_AV ); 

for ( i« (MOV_AV/2 ) +1 ; i< ( f f tdim/2 ) - { MOV_AV/2 ) -1 ; i++ ) { 

magg_ma [ i ] s magg_ma ( i-1 ] ; 

magg_ma[i] -= ((magg[i - (MOV_AV/2) - 

1 J )/( double ) MOV_AV ) ; 

magg_ma[i] += ( (magg[ (MOVJW/2 ) + i] )/ (double )MOV_AV) ; 

} 
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/* within prescribed 'scale zone', calculate final number for 

this angle */ 

angle_vector(angle_int] = 0.0; 

for(i»START_SCALE_ZONE;i<(fftdim/2) - (MOV_AV/2) -l;i++){ 
S mag * magg ( i J / magg_ma { i ] ; 

if( mag > MAGG_THRE S HOLD ){ 

mag -» MAGG_THRESHOLD; 

angle_vector [ angle_int J +« (mag* mag); 

} 

10 , } 

> 

total_angles = angle_ int; 

/* sort out the TOP_CANDIDATES and find which has the best match on 
15 absolute scale */ 

/* chhose the highest angle_vector number for starters */ 
highest = 0.0; 

for ( angle_ int=0 ; angle_int<total_angles ; angle_int++ ) { 
if ( angle_vector [ angle_int ] >highest ) { 
20 highest = angle_vector [angle_int ] ; 

top_candidate = angle_int; 

} 

} 

25 printf ("\n\n tilt from original found = %d (top_candidate/2) - 

45); 

coss a cos(ANGLE_INCREMENT * (double) top_candidate) ; 
sinn » s in ( ANGLE — INCREMENT * (double) top_candidate) ; 
30 /* fill radial vector for this angle */ 

for(i=0;i<(fftdim/2) ;i++) { 

radius - (double) i; 

x » (double) center - radius * coss; 

y = (double) center - radius * sinn; 
35 mag = get_mag(x f y) ; 

x - (double) center + radius * sinn; 

y = (double) center - radius * coss; 

mag +» get_mag(x,y ) ; 

magg[i] = mag; 
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} 

/* create moving average */ 
magg_maIMOV_AV/2] « 0.0; 
for { i»0 ; i<MOV_AV; i++ ) { 
5 maggjna ( MOV_AV/ 2 J +- magg [ i ] ; 

} 

magg_ma(MOV_AV/2] /» ( ( double )MOV_AV ); 

for ( i= (MOV_AV/2 ) +1 ; i< { f f tdim/2 ) - (MOV_AV/2 ) -1; i++) { 

magg_ma ( i J = magg_ma [ i- 1 ] ; 
10 magg_raa{ij -= ((magg[i - (MOV_AV/2) - 1 ] ) / ( double )MOV_AV) ; 

magg_ma[i] += ( (magg( (MOV_AV/2 ) + i] ) / ( double )MOV_AV) ; 

} 

/* now slide the scale and find the highest point */ 
highest * 0.0; 

15 for(scale = SCALE_START; scale < SCALE_STOP; scale+=SCALE_STEP) { 

mag = 0.0; 

f or ( j =0 ; j <number_f r eqs ; j ++ ) { 

radius = scale * freq(j) * (double) fftdim / PI / 2.0; 
if( (int)radius <=* (l+MOV_AV/2) \\ ( int ) ( radiue+1) >* 
20 ((fftdim/2) - (MOVJW/2) - 1) ); 

else { 

frac - radius - (double) ( (int) radius); 
mag +» (1.0-frac)*(magg[ (int)radius ] / 

magg_ma [ ( int ) radius ] ) ; 
25 mag +» frac*(magg[ ( int) (radius+1) ] / 

magg_ma[ (int) (radius+1) ] ) ; 

> 

} 

if (mag > highest) { 
30 highest * mag; 

highest_scale » scale; 

} 

> 

printf ( H \n\nscale found - %lf \n w , 1 .0/highest_scale) ; 

35 

/* now find the exact offset and orietnation, i.e. if it is flipped, 

etc. */ 
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/* display the result at correct rotation and pixel size */ 
sleep (1000) ; 



/* free and clean up */ 
if (DISPLAY_IT) gexit(); 
free(img) ; 
free(disp) ; 
10 free(ar); 

free(ai) ; 

if ( channels 3 3) { 

f ree( img_r) ; 
free( img_b) ; 

15 } 
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